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We present a theory of bilayer two-dimensional electron systems that host a spatially indirect 
exciton condensate when in thermal equilibrium. Equilibrium bilayer exciton condensates (BXCs) 
are expected to form when two nearby semiconductor layers are electrically isolated, and when the 
conduction band of one layer is brought close to degeneracy with the valence band of a nearby 
layer by varying bias or gate voltages. BXCs are characterized by spontaneous inter-layer phase 
coherence and counterflow superfluidity. The bilayer system we consider is composed of two tran¬ 
sition metal dichalcogenide monolayers separated and surrounded by hexagonal boron nitride. We 
use mean-field-theory and a bosonic weakly interacting exciton model to explore the BXC phase 
diagram, and time-dependent mean-field theory to address condensate collective mode spectra and 
quantum fluctuations. We find that a phase transition occurs between states containing one and 
two condensate components as the layer separation and the exciton density are varied, and derive 
simple approximate expressions for the exciton-exciton interaction strength which we show can be 
measured capacitively. 

PACS numbers: 71.35.-y, 73.21.-b 


I. INTRODUCTION 

Recent advances in the study of two-dimensional van 
der Waals material^ have opened up new horizons in 
condensed matter physics by allowing familiar proper¬ 
ties, including those of metals, superconductors, gap¬ 
less semiconductors, semiconductors, and insulators, to 
be combined in new ways simply by designing stacks 
of atomically thick layers. In this article we consider 
condensation of spatially indirect excitons in the case of 
a two-dimensional semiconductor bilayer formed by two 
group-VI transition metal dichalcogenides (TMD) that 
are separated and surrounded by an insulator, for ex¬ 
ample hexagonal boron nitride (hBN). The TMDs are 
in their 2H structure monolayer form. Two-dimensional 
material stacks of this type are promising hosts for exci¬ 
ton condensation, both because they host strongly bound 
excitonsand because of recent progress in realizing 
flexible high quality TMD heterostructures 

In van der Waals heterostructures it is possibld^^^ 
to tune the positions of the Fermi levels in individual 
layers over wide ranges while maintaining overall charge 
neutrality, either by applying a gate voltage between sur¬ 
rounding electrodes or a bias voltage between the semi¬ 
conductor layers. When the indirect band gap between 
the conduction band of one layer and the valence band 
of the other layer is reduced to less than the indirect ex¬ 
citon binding energy, charge will be transferred between 
layers in equilibrium. At low densities, the transferred 
charges form spatially indirect excitons, and these are 
expectecP^Jf^ to form Bose condensates. The bilayer ex¬ 
citon condensate (BXC) state has spontaneous interlayer 
phase cohere nce an d supports dissipationless counterflow 
super cur rent that could enable the design of low- 

dissipation electronic devices 

Exciton condensates in TMD heterostructures are sim¬ 


ilar to atomic spinor Bose-Einstein condensates because 
of the presence of both spin and valley degrees of free¬ 
dom. The spin-valley coupling of conduction band elec¬ 
trons and valence band holes that are specific to TMD 
heterostructure^^ enriches the excitonic physics. In this 
paper we study the interplay between the exciton conden¬ 
sation and spin and valley internal degrees of freedom to 
construct an exciton condensate zero temperature phase 
diagram as a function of effective layer separation d and 
exciton chemical potential /i, or equivalently exciton den¬ 
sity. We demonstrate that there are two distinct conden¬ 
sate phases with different number of condensate flavors, 
as shown in Fig. 

Our paper is organized as follows. In Sec. [llj we ex¬ 
plain how we model the heterostructure, and present the 
mean-field phase diagram implied by Hartree-Fock the¬ 
ory. In Sec. |III[ we derive an effective boson model that 
incorporates exciton-exciton interaction effects and can 
be used to describe excitons in the low density limit. The 
difference between the strengths of the repulsive interac¬ 
tions between excitons with the same internal label and 
between excitons with different internal labels changes 
sign as the layer separation increases. This change drives 
the transition from phase-II, a phase with two conden¬ 
sate flavors present, to phase-I, a phase with only one 
condensate flavor. Both phases spontaneously break the 
symmetry of the model Hamiltonian, and the symme¬ 
try breaking pattern of each phase is analyzed. In this 
section we also explain how capacitance measurements 
can be used to study the exciton phase diagram experi¬ 
mentally and to extract the value of the exciton-exciton 
interaction strength within each phase. In Sec. |IV[ we 
use a time-dependent Hartree-Fock theory to study the 
stability of phase-I against small fluctuations, and to cal¬ 
culate the collective mode spectra of these exciton con¬ 
densates. Finally in Sec. [V| we present a brief summary. 






2 



0 0.2 0.4 ^, ^ 0.6 0.8 1 

d / q,q 


FIG. 1: (Color online) Transferred charge density and ex- 
citon condensate phase as a function of effective layer sep¬ 
aration d and chemical potential fi. The solid orange line 
marks the second-order phase transition from the phase in 
which there is no charge transfer between the bilayers i. e. the 
phase with no excitons present, to the bilayer exciton conden¬ 
sate (BXC) phase. The blue dashed line separates phase-II 
in which two condensate flavors are present from phase-I in 
which the ground state has a single condensate flavor. The 
two red dotted lines are contours at the density values — 
0.01 and 0.1. (For na*B ^0.1 the exciton-condensate ground 
state is expected to be superseded by an electron-hole plasma 
state. See text for a more complete discussion.) The ar¬ 
rows mark the values of d examined in Fig. ia) and [^b), 
and the crosses (:^ indicate the parameter values examined 
in Fig. I^c) and|^d). The left inset is a schematic experi¬ 
mental setup for BXC studies in which the spatially indirect 
gap is tuned by an interlayer bias potential, and the right in¬ 
set illustrates the bilayer band structure in the absence of the 
bias potential 14. 


discuss issues related to experiments, and comment on 
the relationship between our work and previous studies. 


II. MEAN-FIELD PHASE DIAGRAM 


We consider two monolayer TMD semiconductors sep¬ 
arated and surrounded by hBN (Fig. Q. Many of the 
points we make apply with minor modification, how¬ 
ever, to any bilayer two-dimensional semiconductor sys¬ 
tem. Monolayer TMDs are direct-gap semiconductors 
with band extrema located at valleys K and K'. Because 
these TMD layers lack inversion symmetry, spin degener¬ 
acy in the TMD bands is lifted by spin-orbit interactions. 
Because of differences between the orbital character of 
conduction and valence band states,!^ it turns out that 
spin splitting is large at the valence band maxima and 
small at the conduction band minima. As illustrated in 
Fig.[TJ we therefore retain in our theory the two valley- 
degenerate valence bands labeled by v = 1,2, and four 


conduction bands with labels c = 1,2, 3,4 correspond¬ 
ing to spin and valley. We assume that exciton bind¬ 
ing energies and densities are small enough to justify a 
parabolic band approximation for all band extrema. Our 
mean-field ansatz allows up to two types of excitons to 
be present; for examples pairs formed from holes in band 
V = 1 and electrons, selected by spin-splitting, in band 
c = 1, can condense, along with pairs formed from holes 
in band v = 2 and electrons in band c = 2. Although 
the unpaired conduction bands c = 3,4 are only slightly 
higher in energy, this pairing ansatz is fully self-consistent 
at low exciton density, because of the substantial exciton 
binding energy. Our pairing ansatz is also justified by 
an interacting boson model, described in Sec. |ml which 
allows for the most general possible pairing scenario. 

The ansatz leads to the mean-field Hamiltonian: 


Hmf 


k{vc) 









+ E 

k,c=3 

where the prime in the first summation restricts the pair 
index {vc) to (11) and (22) contributions, and a are 
fermionic creation and annihilation operators. The ki¬ 
netic term = h‘^k^{l/{4me) — l/(4m/j,)) accounts for 
the difference between conduction and valence band ef¬ 
fective masses, rUe and ruh^ and are Pauli matrices. 
The dressed energy difference between conduction and 

valence bands, and the coherence induced effective 

k 

interlayer tunneling amplitude, are defined as: 
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( 2 ) 


where (...) is the expectation value in the mean-field 
ground state. 


^(vc) _ J Jvc)2 A (vc)2 

In Eq. m = memh/{me-\-mh) is the reduced mass, 
and A is the area of the system. The paramter jl is: 

= Ane^nd/e, ^ E E AFck '>’ 

k c=l,2 

where /i is the chemical potential for excitons, and n 
is the total charge density transferred between layers. 
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FIG. 2: (Color online) (a)-(b) Energy difference {5E = 
E\\ — E\) per area between phase-I and II states as a func¬ 
tion ji at d/as — 0-1 (a) and d/a^ = 0.5 (b). (c)-(d) Typi¬ 
cal quasiparticle energy bands in phase-II (c), and in phase- 
I (d) with (blue) and without (black) interlayer coherence, 
(d/a^,is (0.1,-2.5) in (c) and (0.5,-1.1) in (d), corre¬ 
sponding to the two crosses(x) in Fig. These results were 
calculated with rUe = ruh- 


Equations (U, (§ and Q form a set of mean-field equa¬ 
tions that can be solved self-consistently. Note that the 
( 11 ) and ( 22 ) pairing channels are coupled through the 
dependence of jl on the total transferred density n. 

The exciton chemical potential can be tuned electri¬ 
cally by applying a bias potential VI between the electri¬ 
cally isolated layers: /i = Vi^ — Eg where Eg is the spa¬ 
tially indirect band gap between the conduction band of 
the electron layer and the valence band of the hole layer. 
The band gap Eg can be adjusted to a conveniently small 
value by choosing tw o-dim ensional materials with favor¬ 
able band alignment 

= 27Te^l{eq) and are the 

Coulomb interaction potentials within and between lay¬ 
ers. The forms of Coulomb potentials are determined 
by solving the Poisson equation for our schematic ex¬ 
perimental setup(Fig. fTl). e = where and ey 

are hBN dielectric constants perpendicular and parallel 
to the z-axis, is the effective dielectric constant due to 
insulator layer (hBN) between electron and hole layers. 
d = D^e^/e\\ , where D is the geometric layer separa¬ 
tion between electron and hole layers, is the effective layer 
separation and slightly larger than 

Below we express lengths and energies in terms of 
the characteristic scales = eh?/{me^)^ and Ry* = 
e^/(2ea^). Typical values for different material combi¬ 
nations are listed in Table HI 

The indirect exciton binding energy E}j determines the 
value for /i at which excitons first appear. When /i < 
—El, no excitons are present. In this state each layer is 
electrically neutral and there is no interlayer coherence. 
Eq. ([2| has nontrivial (n, ^ 0) solutions only for /i > 

—\Ve find two distinct types of BXC phase. In phase- 


TABLE I: Parameter values for different combinations of 
monolayer layer 2H-TMDs. mo is the electron bare mass, 
e = 5 for hBN. Electrons reside in M 0 S 2 , and holes in other 
TMDs. The listed energy gaps Eg apply in the absence of a 
gate voltage. 



me/mo [32] 

m/^/?T^o[32] 

as(A) 

Ry* (meV) 

g.(eV'll^ 

MoS2/MoTe2 

0.47 

0.62 

9.89 

145 

1.1 

MoS2/WSe2 

0.47 

0.36 

12.97 

111 

1.4 

MoS2/WTe2 

0.47 

0.32 

13.89 

104 

0.8 


I, only one type of exciton condenses {e.g. 7 ^ 0 and 

^ = 0). In phase-II, excitons associated with both 
valence bands condense and have equal population (^. 
A^^^^ = A^^^^ 7 ^ 0). Both phases are allowed by Eq. 0. 
We obtain the phase diagram in Fig. by comparing 
the total energy of phase-I and II as a function of (d, ja ). 
Below a critical layer separation dc ~ 0.25a^, phase-II 
always has a lower energy, as illustrated in Fig. [^a). 
Above dc, a transition from phase-I to phase-II occurs 
as the chemical potential /i increases (Fig. |^b)). Typical 
quasiparticle energy bands in phase-II and I are depicted 
in Fig. [^c) and (d), and show that the system is an 
excitonic insulator with a charge gap. 

In our mean-field theory, condensation of one type 
or the other always occurs at T = 0 when excitons 
are present. It is well known however that at high 
electron and hole densities a first-order Mott transi¬ 
tion occur^^^H^ from the gapped exciton condensate 
phase to an ungapped electron-hole plasma state. The 
electron-hole plasma state is preferred energetically be¬ 
cause it can achieve better correlations between like- 
charge particles, reducing the probability that they are 
close together, while maintaining good correlations be¬ 
tween oppositely-charged particles. The density at which 
the Mott transition occurs is most reliably estimated 
via a non-perturbative approaches.!^ No estimate is cur¬ 
rently available for the TMD case, for which the valley 
degeneracy and the small spin-splitting in the conduction 
band will tend to favor plasma states over exciton con¬ 
densate states. Based on existing estimateJ^we can con¬ 
clude that the Mott transition density is below ^0.3 
as d/a^ ^0 and below ^ 0.05 for d/a^ ^ 1 . Cor¬ 
rections to mean-field theory which go in the direction 
of favoring plasma states can be partially captured by 
accounting for screening of the electron-hole interaction 
which becomes stronger as exciton sizes increase and ex¬ 
citons correspondingly become more polarizable. The re¬ 
sults reported here are intended to be reliable only in the 
low exciton density limit. 


III. INTERACTING BOSON MODEL 

To understand the phase diagram more deeply, we em¬ 
ploy a boson Hamiltonian designed to describe weakly- 
interacting excitons in the low density limit. Our strategy 
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FIG. 3: (Color online)(a)Chemical potential /x as a function 
of density n obtained from Hartree-Fock calculations in which 
a single exciton flavor is condensed. The line is a linear ht of 
the numerical data to Eq. (b)Exciton-exciton interaction 
strength extracted from the self-consistent Hartree-Eock equa¬ 
tion solutions (blue dashed line), and exciton-exciton inter¬ 
action strength calculated from the interacting boson model 
(red solid line). The black dashed line is the Hartree contri¬ 
bution qh to the total interaction strength. (c)The blue and 
red lines separate the interlayer (Qx^) and intralayer ex¬ 

change contributions to the total exciton-exciton interaction 
strength. 


to obtain the boson Hamiltonian is to construct a La- 
grangian based on a al wavefunction which parametrizes 
a family of states with electron-hole coherence. The 
Berry phase part of the Lagrangian has the same form as 
that in the field-theory functional integral representation 
of a standard interacting boson model .1^ Appealing to 
this property, we promote variational parameters in the 
wavefunction to bosonic operators. The details of the 
derivation are presented in Appendices and 
The boson Hamiltonian is: 

( 5 ) 

where is a bosonic operator for an exciton with 

a hole in valence band an electron in conduction 
band c, and total momentum Q. Qab is the momen¬ 
tum transfer Qa — Qb- For the TMD system, there 
are 8 possibilities for the composite index {vc). The 
quadratic term in Eq. accounts for exciton kinetic 
energy (M = rrie + 'irih) and chemical potential. The 
quartic terms describe exction-exciton interactions. The 
prime on the quartic term summation enforces momen¬ 
tum conservation Qi + Q 2 = Qs + Q 4 - 

The two types of exciton interaction arise from the 
fermionic Hartree and exchange interactions respectively. 
In the exchange interaction, two excitons swap con¬ 
stituent electrons or holes. Analytic expressions for the 
coupling strength gni^ and are given in Ap¬ 


pendix 

We focus here on their zero-momentum limits gn = 
and gx = which are more easily inter¬ 

preted and capture much of the exciton-exciton interac¬ 
tion physics. For the case in which the exciton conden¬ 
sate is populated by a single flavor we find that for low 
exciton densities 


11 = -Eb + gn. 


(6) 


where g = gn E gx the total exciton-exciton in¬ 
teraction, as expected from the mean-field theory for 
weakly interacting bosons. This behavior is illustrated 
in Fig. [^a). We have verified that the interaction pa¬ 
rameter obtained by examining the dependence of g on 
n in the fermion mean-field theory agrees with the an¬ 
alytic expression in App. as illustrated in Fig. |^b) 
which plots ^ as a function of layer separation d. We find 
that gn = EKC^dje and that gx = g^x ~ 9^x where g^^ 
and g^x both positive and originate from inter and 
intra layer fermionic exchange interactions respectively. 
The binding energy of isolated excitons is due micro¬ 
scopically to attractive inter layer exchange interactions. 
When excitons overlap and interact with each other, co¬ 
herence between layers is reduced weakening inter layer 
exchange, but strengthening intra layer exchange. This 
explains the signs of the two contributions to gx- The 
overall sign of gx is positive at d = 0 because the loss 
of interlayer exchange energy when excitons overlap is 
greater than the gain in intralayer exchange energy. In 
Fig. [^c), we show that gx becomes negative beyond a 
critical layer separation dc ^ 0.25a^. It turns out that 

although both g^^ 9^x increase with layer separation 

d, the rate of increase of g^^ i^ smaller than for g^^ • The 
difference in behavior can be traced to the exponential 
decrease in the momentum space inter-layer Coulomb in¬ 
teraction with layer separation d as shown in Eq. (A19) 
and (A20). 


To find the ground state in the realistic multi-flavor 
case, we assume that all excitons condense into (5 = 0 
states and introduce the following matrix: 


jr _ A-®(11)) (-^(12)) (-^(13)) (-6(14))) 

\/4 ((- 8 ( 21 )) {B(23)) {B(2i))j 

Neglecting the small spin-orbit splitting of conduction 
band states, the total energy per area can be written in 
a compact form, 

= -{Eb + M)Trr + ^-{TvTf + ^TrT^, (8) 

where T = In Eq. ^ TrT is the total density of 

excitons,summed over all flavors,and TrT^ — (TrT)^ mea¬ 
sures the flavor polarization of the exciton condensate. 
This energy functional is invariant under the following 
transformation: 


T^E'^U2BUI. 


( 9 ) 
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Here captures the U(l) symmetry which originates 
from separate charge conservation in the individual lay¬ 
ers. IA 2 and Z //4 are respectively 2 x 2 and 4x4 spe¬ 
cial unitary matrices, which capture the SU( 2 ) symme¬ 
try of the valence bands and the SU(4) symmetry present 
in the conduction bands when their spin-splitting is ne¬ 
glected. The overall symmetry group of the system is 
U(l) xSU( 2 ) xSU(4). When the conduction band spin- 
orbit splitting is included, the higher energy conduction 
band states in each valley are not occupied and the sym¬ 
metry group is reduced to U(l) xSU( 2 ) xSU( 2 ), corre¬ 
sponding to separate charge conservation and rotations 
in both conduction and valence band valley spaces. 

T acquires a nonzero value in the ground state only 
if /i > —£’ 5 , . By minimizing the energy functional, we 
verify that the sign of gx determines the position of a 
phase boundary between two different classes of exciton 
condensate which we refer to as phase-I and II. When 
< 0, phase-I is energetically favorable and a repre¬ 
sentative realization of the ground state is. 




fl 0 0 o\ 

yo 0 0 oy ’ 


( 10 ) 


where ni = ^gx) is the exciton density. 

is invariant under the transformation: 


0 




(11) 


where V 3 is a 3 x 3 unitary matrix. Therefore, phase-I 
spontaneously breaks the U(l) xSU( 2 ) xSU(4) symmetry 
down to U(l)xU(3) symmetry. 

When gx > ^ phase-II is realized. Energy minimiza¬ 
tion shows that a representative realization of the ground 
state in phase-II is. 


= ,AW 2 (j J ° 2 ) . ( 12 ) 

where nn = (/i + Ei))/{gH + ^x/ 2 ) is the total exciton 
density in phase-II. Eu is invariant under the transfor¬ 
mation: 

where V 2 is a 2 x 2 unitary matrix. Phase-II sponta¬ 
neously breaks the U(l) xSU( 2 ) xSU(4) symmetry down 
to SU( 2 ) xU( 2 ) symmetry. A similar analysis can be ap¬ 
plied to identify the symmetry breaking pattern when 
the spin splitting of the conduction bands is consid¬ 
ered. In phase-I, an application of an infinitesimal ex¬ 
ternal Zeeman field lifts both conduction and valence 
band valley degeneracies, and selects a unique conden¬ 
sate ground state with a finite spin-polarization, as illus¬ 
trated in Fig.[^ Phase-I therefore satisfies the definition 
of a ferromagnet, defined as a system with a finite spin- 
polarization in an infinitesimal Zeeman field, but has a 


(a) (tt) 
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FIG. 4: (Color online) (a) Four representative degenerate 
ground states in phase-I. For illustration purpose, only two 
conduction bands are shown. (The other two are slightly 
higher in energy because of spin-orbit splitting and do not 
participate in the ground state manifold.) The dashed black 
oval highlights the two bands with spontaneous phase co¬ 
herence in phase-I. Charge is transferred from the valence 
band partner to the conduction band partner. In the upper 
panels (tt) £ind (ii), coherence is established between like 
spins and the ground state is not spin-polarized. In the lower 
panels (ti) and (it)? charge is transferred between opposite 
spins and the ground state is spin-polarized, with a finite 
spin-polarization that is proportional to the transfered charge 
density. Because spin and valley are locked by spin-orbit cou¬ 
pling, spin-polarized states are stabilized by an infinitesimal 
Zeeman field, (b) In the presence of an infinitesimal Zee- 
man field that favors spin up, the spin-polarized state (^t) is 
selected as the unique ground state in phase-I. (c) Two rep¬ 
resentative degenerate ground states in phase-II. Both states 
have zero spin-polarization, and remain degenerate in an in¬ 
finitesimal Zeeman field as shown in (d). 


distinct set of broken symmetries compared to the usual 
spin rotational symmetry breaking. 

Based on this mean-field calculation, we conclude that 
although the system has 8 types of excitons in total, only 
one or two flavors condense in the ground state. The 
number of condensed flavors is in general limited by the 
number of distinct valence or conduction bands, which 
ever is smaller in number. Although the boson model 
correctly captures the phase transition position as a func¬ 
tion of d, it is important to emphasize that it is valid 
only in the low exciton density limit. For this reason, it 
fails to accurately predict the g dependence of the phase 
boundary. In addition it fails to capture the tendency 
toward weaker electron-hole pairing at high exciton den¬ 
sities, which eventually leads to an electron-hole quantum 
liquid state with no interlayer coherence. 

The relationship between the exciton density n, the 
exciton chemical potential g and the coupling strength 
g makes it possible to extract the value of g from ca¬ 
pacitance measurement. The differential capacitance per 
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area for the heterostructure is: 

Ofi 


(14) 


The Hartree coupling strength qh can be identified as 
the inverse of the geometric capacitance: 


Cgeo = (? I Qh = e/{47rd). 


(15) 


Therefore, capacitance measurement provides a simple 
way to determine the value of gx in the low-exciton den¬ 
sity limit: 


^ phase-II 

The sign of gx helps to distinguish phase-I and IL 


S'x < 0, phase-I 


(16) 


IV. FLUCTUATIONS AND STABILITY 

The bilayer exciton condensate is a state with spon¬ 
taneously broken continuous symmetries, and therefore 
hosts low-energy collective fluctuations. Theoretical 
studies of fluctuation properties are of interest in part 
because they can reveal mean-field stat^^ instabilities. 
The collective modes can be studied using the interact¬ 
ing boson model, which is described in detail in Ap¬ 
pendix The interacting boson model admits analytic 
solutions for collective modes associated with exciton 
density, phase and flavor fluctuations in both phase-I and 
II. However, it is valid only in the low exciton density 
limit. Here, we study another approach that can be ap¬ 
plied to any exciton density. This approach is based on 
the following variational wave function which captures 
exciton density and phase fluctuations in phase-I: 


i*>=n [%+E J 

k Q 

vt 


(17) 


where \XC) is the phase-I ground state. 7 I ^ and 7 I ^ are 

respectively quasiparticle creation operators for occupied 
and empty quasiparticle states in \XC) associated with 
the ground state condensate and are defined as follows: 




7I =uto} ^ 

'k,0 k ck 

7I =Vra^^ 

k,l k k 


(18) 


where and are parameters determined by self- 
consistent Hatree-Fock equations 


= l7l-4V%), n = \Ui^+^k/Ej:)- (19) 


Zj^ is a normalization factor. 






- eh continuiin 
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' phase 


0.0 0.5 1-On^* 1-5 2.0 





FIG. 5: (Color online) (a)-(b) Eigenvalues of for d/a% = 
0.5 as a function of momentum Q at na*B = 0.008 (a) and 
= 0.08 (b). Red solid and blue dashed lines connect 
the lowest eigenvalues of the density fluctuation matrix 
and the phase fluctuation matrix }C^~^ respectively. Small 
red points label higher eigenvalues of The gray lines 

mark the lower edge of the quasiparticle electron-hole con- 
tinua [Min(E^ + (c)-(d) Collective mode spectra at 

— 0.008 (c) and na*B = 0.08 (d) for = 0.5. Larger 
points are used for the lowest energy excitations. In (c), the 
red and black lines are respectively a linear fit at small Q and 
a quadratic fit at large Q to the collective mode energy. For 
the results presented here, we assumed rUe = ruh- 


and z^{Q) are complex parameters. 

To study fluctuation dynamics, we construct the La- 
grangian: 


c = B 


( 21 ) 


where B is the harmonic Berry phase, and is the 

harmonic energy variatiorP^ 


5^(2) = {S-^^p{Q)zl{Q)zp{Q) 

Q,k,p 

+ h,-p(Q) 7 (e)^p(-Q) + ^|(Q) 7 (-Q)]}- 


( 22 ) 


Explicit forms for the matrices £ and T are give n in 
App. 0 To decouple ZQ contributions in Eq. ( |22| , we 
perform a change of variables, defining 


^kiQ) ~~ ^ ’WkiQ) 


(23) 

(24) 


Note that x^{Q) and y^{Q) are complex numbers, and 
that there is a redundancy. 


x_i:{-Q) =x*^(Q), 

V-ki-Q) = vliQ)- 


(25) 


( 20 ) 


In terms of the x and y fields, the Berry phase B and 
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energy variation are, 

B = i ^ (y*dtx + ydtx* - x*dty - xdty' 




Q,k 


(5£'(2) - + 2/|^^j2/p)(Q)- 


(26) 


Q k,p 


The kernel matrices (Q) = -p) (^) 

and symmetric. The x and y fields in can be identi¬ 

fied with exciton density and phase respectively, and the 
Berry phase contribution to the action captures the con¬ 
jugate relationship between these fluctuation variables. 

Stability of the mean-field ground states against small 
fluctuations requires that the matrices are nonnega¬ 
tive. We have verified that this condition is satisfied out 
to large d by explicit numerical calculations like those 
summarized in Fig. [^a) and (b). At Q = 0, the matrix 
JC^~^ always has a zero-energy eigenvalue since, 

p P 


which follows from the fact that ground state energy is 
independent of global interlayer phase. 

For low exciton density n (Fig. [^a)), the lowest eigen¬ 
values of and JC^~^ have similar behavior and are 
separated from the continuum. This is expected since 
and IC^~^ are identical in the limit n ^ 0. Fig. [^b) 
demonstrates that the lowest eigenvalues of are close 
to the particle-hole continuum when the exciton den¬ 
sity becomes large; the interacting boson model discussed 
above fails qualitatively in this limit. 

The Euler-Lagrange equation for the Lagrangian in 
Eq. (21) gives rise to the equation of motion. 


dtXf:{Q) = -{lC^^^.yp){Q), 
dty^,iQ) = +{lC^^Up)iQ). 


(28) 


which leads to 


dhdQ) = - [ («)• 


(29) 


It follows that the energy of the collective mode is given 
by the square root of the lowest eigenvalues of the matrix 
product which is plotted in Eig. [^c) and (d). 

The lowest energy collective mode is the gapless Gold- 
stone mode of the exciton condensate. Eor low exciton 
density n (Eig. |^c)), the Goldstone mode has linear dis¬ 
persion at small Q, becoming quadratic at large Q. This 
agrees with the Goldstone mode behavior predicted by 
the weakly interacting boson model(Eq. (§)• Eor large n 
(Eig.j^d)), the Goldstone mode deviates from quadratic 
behavior at large Q. The failure of the weakly interact¬ 
ing boson model in the high density limit originates from 
the internal structure of the excitons. When the typical 
distance between excitons is comparable to exciton size, 
excitations must be described in terms of th e under lying 
conduction and valence band fermion states.^ ^^ * ^^ * 


V. SUMMARY AND DISCUSSION 

By combining Hartree-Eock theory and an interacting 
boson model, we have shown that spatially indirect ex¬ 
citon condensates in group-VI TMD bilayers have two 
distinct phases. We have also studied the dynamics of 
exciton condensate density and phase fluctuations and 
calculated the associated collective mode spectra. 

The topic of exciton condensation in semiconductors 
has a long history and our work is related to some earlier 
studies. Eor example, Berman et al^ studied exciton 
condensation in bilayers formed from gapped graphene, 
although the possibility of two distinct condensate phases 
was not considered. The phase transition between the 
two condensate phase s as a function of layer separation 
was studied previousl}BI^ for the case of quantum well 
bilayer excitons, and further explored in a very recent 
publication.!^ The TMD layers considered in this pa¬ 
per are distinguished from semiconductor quantum well 
systems by exciton binding energies that are an order 
of magnitude larger, and by spin-valley coupling which 
leads to two-fold degenerate valence bands and approx¬ 
imately four-fold degenerate conduction bands. Com¬ 
pared to Refs|45] and |46l we used a completely differ¬ 
ent approach to derive an interacting boson model. Our 
approach is physically transparent, and is based on a 
variatonal wavefunctions deflned by parameters whose 
quantum fluctuations are characterized by using a La¬ 
grangian formalism. The bosonic nature of excitons is 
automatically taken into account in the Lagrangian, and 
there is no need to calculate combinatorial factors aris¬ 
ing from the indistinguishability of bosonic particles. Our 
approach provides a simple yet systematic way to model 
the exction-exciton interaction. We have also discussed a 
fermionic Hartree-Eock approach from which the exciton- 
exciton interaction strengths can be extracted with sim¬ 
ilar results, and proposed that the interaction strengths 
can be experimentally determined by performing capac¬ 
itance measurement. 

Because the hBN dielectric barrier in the systems of 
interest, must be thick enough to make interlayer tun¬ 
neling weak, Eig. implies that phase-I with a single 
condensate flavor is more likely to be realized in experi¬ 
ment than phase-II. Phase-I breaks the invariance of the 
system Hamiltonian under separate valley rotations in 
conduction and valence bands, and is ferromagnetic in 
the sense that inflnitesimal Zeeman coupling leads to a 
spin-polarization that is proportional to the exciton den¬ 
sity. 

In spit of their large gaps, band edge states in TMDs 
have relatively large Berry curvatures^. In monolayer 
TMDs momentum space Berry curvatures leacP to un¬ 
usual exictonic spectra in which hydrogenic degeneracies 
are lifted. Although band Berry curvatures should be less 
important in spatially indirect exciton systems because 
weaker binding implies that the exciton states are formed 
within a smaller region of momentum space. In terms of 
its influence on quasiparticle bands, exciton condensation 









has the effect of preventing gaps between conduction and 
valence band states from closing. Since the host semicon¬ 
ductor materials are topologically trivial, and since tran¬ 
sitions between trivial and non-trivial states can occur 
continuously only when the quasiparticle charged excita¬ 
tion energy vanishes, we argue that exciton condensation 
will not result in interaction-driven topologically nontriv¬ 
ial states in our system. 

The critical temperature of spatially indirect exciton 
condensate is the Berezinskii-Kosterlitz-Thouless transi¬ 
tion temperature, given at low exciton densities by the 
weakly-interacting boson expression 

n TYi 

^ 2.6(na^^)^Ry*, (30) 

where M is the electron-hole pair total mass and m 
is the reduced mass. In the low exciton densit y limit 
scales linearly with exciton density.^^^EZl For the 
MoS 2 /hBN/MoTe 2 heterostructure and exciton densities 
n in the range, 0.01 and Tbkt is about 

lOK. The maximum possible transition temperature is 
closely related to the critical density at which the Mott 
transition to an electron-hole plasma occurs, and this in¬ 
creases as d/a^ decreases. Usin g t he variational Monte 
Carlo estimate of DePalo et the critical value of 
^ 0.3 as d/a^ ^ 0. From this we conclude that 
^bFbkt cannot exceed around 300K. Adjustment of ex¬ 
citon density by external bias voltage can be employed 
to search for the highest transition temperature and to 
study the Mott transition to an ungapped electron-hole 
plasma that is expected at high exciton densities. The 
most interesting regime is likely to be the case of very 
small layer separations of which current leakage driven 
through the tunnel barrier by an interlayer bias potential 
might be appreciable, requiring the bilayer to be treated 
as a non-equilibrium system. 

The exciton condensate should be experimentally re¬ 
alizable in TMD bilayers provided that samples with 
sufficiently weak disorder can be achieved. The photo¬ 
luminescence line width W of an individual monolayer 
TMD is a particularly useful characterization of sample 
quality for this purpose. W is currently dominated® by 
the position-dependence of exciton energies. Therefore, 
the narrower the line width VF, the weaker the disor¬ 
der, and the better the sample quality. We expect bi¬ 
layer exciton condensation to occur only in samples in 
which W < A^bFbkt, since the excitons will otherwise 
simply localize near positions where they have minimum 
energy. Note that the inhomogeneous broadening W of 
spatially indirect excitons will not be experimentally ac¬ 
cessible since the corresponding transitions are optically 
inactive when the interlayer tunneling is negligible, but 
that it should be similar to the broadening of the readily 
measurable direct exciton energies. It should therefore 
be possible to judge on the basis of optical characteri¬ 
zation when samples have achieved sufficient quality to 
study spatially indirect exciton condensate physics. 
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Appendix A: Interacting boson model for excitons 
in the low density limit 


In the low density limit, excitons can be approximated 
as interacting bosons. We take a BCS like variational 
wave function to describe excitons. 


1^) = ^exp(n1’)|vac), 

U’*’ = Ayca^gy, 
v,c 


(Al) 


where C and V respectively denote a conduction and 
valence band state, and include internal indices such as 
spin and valley and also momentum label. |vac) is the 
vacuum state defined by ay|vac) = ad vac) = 0 . 
operator creates particle-hole excitations on top of the 
vacuum. A/' is a normalization factor so that (T|T) = 1 . 
Aye is a set of complex variational parameters, which 
are small when the exciton density is low. 

The density matrix with respect to |T) is pa^ = 
(T|aJ^a/ 3 |T), where a and /3 can be conduction or valence 
states. We expand the density matrix to fourth order in 
Aye, 


pvv^ - 

^ ^yy' + (—AA^ + AA^AA^)yy/, 


PCC' - 

^ (A’*'A — A’*'AA’*'A)c'C', 

(A2) 

pvc - 

^ (A — AA^A)yc', 



where A is understood to be a matrix and A^ is its Her- 
mitian conjugate. 

We introduce another matrix A so that pec' has a 
quadratic form without fourth-order correction. 


A = A + -AAU. 


(A3) 


-(a)- 


Iviq) 


iViq) 


>Mq) 


HCT)- 

(a) (b) (c) 


FIG. 6: Feynman diagrams for the Hartree exciton-exciton 
interaction processes. A double line with arrow represents 
an exciton state, and a single solid (dashed) line with ar¬ 
row depicts an electron (hole) in the exciton. Wavy lines are 
interaction or U{^. (a) and (b) are the intralayer con¬ 

tributions, while (c) is the interlayer contribution, (a ), (b) 
and (c) correspond to the three ter ms in (Eq. (A9)) and 
also the three terms in (Eq. (A17)). 
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Expanding p up to fourth order of A, we have that 


Pvv - 

» dvv - (AA'^)vy/, 


Pcc - 

a (A^A)cc', 

(A4) 

Pvc - 

a (A - Aa^A)vc. 


The number of excitons is {Nex) = '^c Pcc ~ TrA’^'A. 
Therefore, we verify that A acts as the small parameter 


in the limit of low (A^ecc)- 

An important property of the density matrix is that 

y - p = 0(A5), (A5) 

which indicates that |^) can be approximated as a Slater 
determinant up to fourth order in AP^ 

1^) parametrizes a family of states with electron-hole 
coherence, and also represents low-energy states in the 
low-exicton density limit. We choose to construct an 
effective interacting boson model using this variational 
wavefunction approach rather than a commonly used 
auxiliary field approach because of the necessity of consis¬ 
tently accounting for both exchange and Hartree mean- 
fields in spatially-indirect exciton systems. (See addi¬ 
tional discussion below.) To study low-energy dynamics, 
we construct a Lagrangian based on |T) 

C = {^idt =B-n, (A6) 

and again expand everything to A^. This Lagrangian 
provides an effective field theory for excitons. The Berry 
phase has the following form, 

B = « ^(lV[AtatA] - IV[(atAt)A]), (A7) 


Um A >. \u{q) 

(a) (b) 



FIG. 7: Feynman diagrams for the exchange exciton-exciton 
interaction processes. The convention is the same as in Fig.[^ 
(a) and (b) are the interlayer contr ibut ions, which correspond 
to the firs t two terms in 1-i^^ (Eq. (A9)) and the two terms in 
(Eq. (A19)). (c) and (d) are the intralayer contribu tions , 
which correspond to the lase t wo te rms in (Eq. (A9)) 

and the two terms in (Eq. (A20)). 


which does not have fourth order corrections. 

To calculate the energy functional 77, we take advan¬ 
tage of the Slater determinant approximatioiPSl to |T) 
(Eq. ( |A5| )) and obtain that 

n = {^H\^) , (as) 

where is quadratic in A, and T-L^px i® quartic in A 
with subscript H and X representing Hartree and ex¬ 
change contributions. The explicit forms are below. 

77^^^ = {^c - ^v)^^cv-^yc 

- Wv^CiC2V2^CiV2^yiC2^ 

'^H = 2^'^1*^2C3C4(A^A)ciC4(A^A)c'2C3 

+ 2^^i^2V'3V'4(AA^)viVi(AA^)v'2V3 

- brviCiC2V'2(AA^)viVb(A^A)ciC2) 

'^x = 

+ 2 ^^iC'iC'2 V2 (A^AA^)ciV2AvjC'2 

- 2^^1020304 (A^A)c'iC3(A^A)c'2C4 

- 2^V"iV2V3V4(AA'^)viV3(AA'^)v2V4. 

Here sc and Sy are conduction and valence state energy 
including self-energy effects. The interaction kernel W 
has the form 


^(ni fci ) (722^2 )(n3fe3) ( 724 ^ 4 ) 

= ^<^721724 ^n2723%i+fc2,fc3+fc4^’^l’^2(^l “ 


(AlO) 


where the momentum dependence is now explicit, and 
n denotes internal indices. A is the area of the system. 
Wnin 2 {^ is the intralayer interaction V{^ if both ni 
and 77-2 represent conduction or valence bands, and the 
interlayer interaction U{q) otherwise. 

We now write 77^^^ in a more concrete form 


77 ( 2 ) 


{v ,k){c,k+Q) 


V 2me 


‘2'mh 



--U{k-k') 


A 


{v,k'){c,k'+Q)' 


(All) 


Here v and c denote different valence and conduction 
bands. We approximate sc and sy by parabolic bands, 
and assume different valence (conduction) bands have the 
same hole (electron) mass rrih {Trie). la the case of TMDs, 
these are reasonable approximations, and v and c respec¬ 
tively take two and four different values. 

77^^^ can be reduced into a diagonal from by doing the 
following decomposition. 


^iv,k)ic,k+Q) ~ +(^ 12 ) 
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where is a complex field that depends on momen¬ 
tum Q but not on k. where the total mass 

M = rUe + rrih- For notation convenience, we also intro¬ 
duce Xe = rrielM = 1 — Xh- f{k) is the Is wavefunction 
for a single exciton, 


. 2m 


^U{k-k')\f{k') = -Ef,fik), 


(A13) 


where the reduced mass m = mem^/M^ and Eij is the 
binding energy for Is state. The normalization condition 
is that 




(A14) 


Here we have chosen f(k) to be real. In Eq. (A12), /(A? + 
XhQ) is the wavefunction for an exciton with center-of- 
mass momentum Q. 

By substituting Eq. (A12) into Eq. (A7) and (A9), we 
obtain that 


S: 




(vc)Q^^^ (vc)Q ■®(wc)Q^*-®(„c)q) 


(A15) 


«=s; (w - 

I ' 

S {vc)Qi^{v' c')Q 2 ^{.Vc')Q3^{,vc)Q4. 

+5'x(Qi3, {^yc)Qi^ {v' c')Q2^ {v' c)Q3^ {vc')Qa) ' 

(Alb) 

The Berry phase B has the same form as that in the 
field-theory functional integral representation of a stan¬ 
dard interacting boson model, which suggests that the 
Lagrangian £ = S — is a functional field integral rep¬ 
resentation of a boson modeP^. By replacing the complex 
numbers with bosonic creation and annihilation 

operators in the energy functional H, we arrive 

at the boson model <§ in the main text. 

gni^ is derived from and has the following ana¬ 
lytic expression 

9h{^ = V{^[F{xe^^ + F{xh^^] 

- 2U{^F{x,^F{xh^, 

where 

k 

Eor zero-momentum transfer, = 2(E(0) — U{Q)) = 

Ane^d/e. 

Q14) is derived from and can be further 

decomposed into two parts gx = — Qx^- Qx^ arises 


from the loss of interlayer exchange energy as more ex- 
citons condense, while —Qx^ comes from the gain of in¬ 
tralayer exchange energy as electron or hole density in¬ 
creases. The explicit forms of are below. 

5'x^(Qi3, Q14) 

=j^E^ik-k')fik')f{k + Qe) 

X [fik + Qh)f{k + Qe + Qh) 

+ f{k-Qh)f{k + Qe-Qh)\ (A19) 

k 

X [/(^ + Qh)f{k + Qe + Qh) 

E f{k — Qh)f{k + Qe — Qh)] 5 


k,k' 

^ ^ ^ ^ \ ^ ^ ^ ^ A20 

V{k -k'E Qe)f{k)f{k')f{k + Qh)f{k' + Qh)E 

Vik -k'E Qh)f{k)f{k')f{k + Qe)f{k' + 4 )}, 


whe re Qp = XqQi^ and Qh = XhQiA- According to 
Eq. (A19), can also be interpreted as the increase of 
kinetic energy due to the increase of the exciton density. 
The diagrammatic representation of different processes 
for the exciton-exciton interaction is shown in Eig. 
and The momentum dependence of the interaction 
strength gH,x is illustrated in Eig. Both gn and gx 
have a strong momentum dependence, which indicates 
that extion-exciton interaction are long-ranged instead 
of short-ranged. 

The problem of exciton-exction interaction has 
been st udied u sing many different methods in the 
hteratur (j^^F6 | 49t t^^ Here we used an alternative ap¬ 
proach based on a variational wave function combined 
with the Lagrangian formalism. The functional field in¬ 
tegral representation of the boson model in Eq. (§ is 
exactly gi ven b y the Lagrangian C with Berr y ph ase 


ry pUase 

|ai^. 


B in Eq. (A15) and energy function 1-L in Eq. 

This property establishes the equivalence between the 
Lagrangian and the boson model. In fact, all results 
obtained using the boson model can be equivalently de¬ 
rived from the Lagrangian. The phase transition between 
phase-I and H can be determined by the minimization of 
the energy functional ( A16| ) with the ansatz that (H*, B) 
is spatially uniform and time independent. The collective 
mode studied in Appendix 0 can be obtained using the 
Euler-Lagrange equation of the Lagrangian. 

There are other approaches in constructing an effective 
theory of bosonic excitation in an interacting fermion sys¬ 
tem. One example is the standard Hubbard-Stratonovich 
transformation. There is however a certain arbitrari¬ 
ness in decomposing electron-electron interactions into 
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01 23 4^ 5678 

qaB 


FIG. 8: gx{q,0) and gx{q,^ as a function of mo¬ 

mentum. d/a^ is 0.1 in (a) and 0.5 in (b). rUe = rrih is 
assumed. 


Hartree or Fock channels in the Hubbard-Stratonovich 
schem^^. The HS approach is not appropriate here be¬ 
cause a proper description of the SIEXC requires that 
Hartree and Fock interactions to be treated on the same 
footing. 


Appendix B: Interacting boson model for excitons 
with zero center-of-mass momentum 


We consider another variational wave function with all 
excitons condense into zero center-of-mass momentum 
state. 

l^-o) = n E ^ 

[x(fc)aG + E hv2c)alf:] |0), 


where vi = 1 and V 2 = 2, representing two valence 
bands, and c represents different co nduct ion bands. f{k) 
is the exciton wave function in Eq. (A13). are com¬ 
plex parameters, which are independent of momentum 
k. In |0), both valence and conduction bands are empty. 


~ Facotrs si^ 2 {k) and x(/c) are de¬ 
termined by normalization and orthogonality conditions. 
By normalization conditions, we have that 


To ensure that at each momentum k the two occupied 
states are orthogonal to each other, we require that 

siik)x{k) + f{kfJ2^Uc)hv,c) = 0. (B3) 

c 

Similar to the procedure in App. we expand the 
energy functional (To|i^|To) to fourth order in 
and then replace complex numbers 

erators which gives rise to the same 

boson model in Eq. ( |A16| ) except that all momenta are 
restricted to be zero. 


Appendix C: Collective modes in the interacting 
boson model 


The collective modes can be calculated analytically us¬ 
ing the interacting boson model. The strategy is to shift 
a bosonic operator by its mean-field value. 




{vc)Q 


— (^Gc)q)^Q, 0 + ^Gc)Q’ 


(Cl) 


where is also a bosonic operator that describes 

fluctuations around the mean field state. The bosonic 
Hamiltonian in Eq. (§ is then expanded to second or¬ 
der in The first order terms vanish as the en¬ 

ergy functional in Eq. is minimized in the mean-field 
state. The quadratic terms can be diagonalized using the 
Bogoliubov transformation, giving rise to the collective 
mode spectra. We present the dispersion and degeneracy 
of collective modes in phase-I and H without giving the 
details of the derivation. 

In phase-I, there are five gapless modes, and three 
gapped modes. Given the mean-field value in Eq. (10), 
{b^,b)(^yc) with different composite index (vc) are decou¬ 
pled. The (vc) = (11) mode is gapless, which is the 
Goldstone mode due to the spontaneously broken U(l) 
symmetry, i.e. the separate charge conservation within 
each layer. Its dispersion is: 


^(11) 




L 2M 


{2gH + ^+)(Q)^i 


\ L 2M 


■g-{Q)ni 


g ±{ Q ) = 5'x(Q, 0) +5'x(0,Q) -^x(0,0) ± gx { Q ’, Q ), 


(C2) 

where g-{Q) has a dependence at small Q. Therefore, 
the (11) mode has a linear dispersion in Q ^ 0 limit. 

The (vc) = (12), (13) and (14) modes are degenerate, 
with a gapless dispersion: 


<^(12) 


2M 


[9x{0,Q) -fl®(0,0)]ni, 


(C3) 


si{k) = ^1- f{k)^J2hvic)\^’ 

S2{k) = ^l-|x(fc)|2-/(fc)2^|6(„,e)P. 


(B2) 


which has a quadratic Q dependence in Q ^ 0 limit. 
Similaryly, the (vc) = (21) mode is also gapless and 
quadratic at small Q: 

W( 21 ) = - 1 ^ + [fe(Q,0) -fla;(0,0)]ni. (C4) 
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The {vc) = (22), (23) and (24) modes are degenerate 
and gapped: 


W( 22 ) = - QxiQ, 0)ni, (C5) 

which shows that phase-I is stable against small fluctua¬ 
tions provided that ^a;(0,0) is negative. 

In phase-II, the mean-field values are given in Eq. 
and all eight collective modes are gapless. The {vc) = 
(11) and (22) fluctuations are coupled, and give rise to 
two non-degenerate gappless modes: 


<^( 11 ),( 22 ) 


(A 


■^( 11 ),( 22 ) 


i 


(C6) 




decoupled, and have degenerate collective modes: 

^(13) = 2^ [9x{0,Q) -^a:(0,0)]^, (C7) 

which has a quadratic dispersion at small Q. 


TABLE II: Classification of collective modes in phase-I and 
II. iVgapped is the number of gapped collective modes. Ni and 
N 2 are respectively the number of gapless collective modes 
with linear and quadratic dispersion. Nbsg is the number of 
the broken symmetry generators. 



Agapped 

Ni 

N 2 

Absg 

phase-I 

3 

1 

4 

9 

phase-II 

0 

4 

4 

12 


In Table [TTj we list Ni and N 2 , respectively the num¬ 
ber of gapless collective modes with linear and quadratic 
dispersion, and A^bsg? the number of the broken sym¬ 
metry generators for each phase. There is a relationship 
among these three numbers in both phase-I and II, 


both of which have a linear dispersion at small Q. 

The (vc) = (12) and (21) fluctuations are also coupled, 
and lead to two modes that are degenerate with ^^ 2 ) • 
The (vc) = (13), (14), (23) and (24) fluctuations are 


A^i -h 2N2 = A^bsg- (C 8) 

This relationship is typical for broken symmetry states 
in system without Lorentz invarianc^^. 


Appendix D: explicit expressions for <^^^(Q) and r^^(Q) 


Below are explicit expressions for £j^ ^{Q) and r^^((3), which appear in the energy variation in Eq. (22). 


k,p 


^k,p(^^ ~ [^(Q) ~ E Q — p)] 

EjU{kEQ W+QVq)’ 


(Dl) 


where and are defined in Eq. (19), and V{Q) and U{Q) are respectively intralayer and interlayer Coulomb 
interactions. ^_ 
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